// TODO: NOTE: The calculation of the bivariate normal distribution should go to functions++.

// TODO: NOTE: The calculation of CDF of distributions, etc. should follow more closely Boost [a general function cdf() that takes a distribution, etc.].

// TODO: NOTE: A better approach to calculating the CDF should be followed.
//
// TODO: NOTE: ... See the references in (for example):
//
//     http://www.itl.nist.gov/div898/software/dataplot/refman2/ch8/bvncdf.pdf


// STL
#include <cmath>                               // std::sqrt(), ::exp()
#include <functional>                          // std::function<>, ::bind(), ::cref()

// Boost
#include <boost/math/distributions/normal.hpp> // boost::math::cdf(), ::normal()

// jScience
#include "jScience/physics/consts.hpp"         // PI

// functions++
#include "functionsxx/integration.hpp"         // functionsxx::integration::quadrature_GaussLegendre()


static double CDF_standardized_bivariate_normal_distribution_outer(
                                                                   const double x,
                                                                   // -----
                                                                   const double rho,
                                                                   // -----
                                                                   const double inf,
                                                                   const double b2,
                                                                   // -----
                                                                   const int    n
                                                                   );


static double CDF_standardized_bivariate_normal_distribution(
                                                             const double x,
                                                             const double y,
                                                             // -----
                                                             const double rho
                                                             );


//
// DESC: Calculate the area under a binormal-LR ROC curve.
//
// INPUT:
//
//     double       d_a : Metz and Pan binormal-LR parameter
//     double       c   : Metz and Pan binormal-LR parameter
//
// OUTPUT:
//
//     double       AUC :
//
// -----
//
// NOTE: See:
//
//    C. E. Metz and X. C. Pan, ``"Proper" binormal ROC curves: theory and maximum-likelihood estimation,'' J. Math. Psych. 43, 1-33 (1999).
//
//    S. L. Hillis, ``Equivalence of binormal likelihood-ratio and bi-chi-squared ROC curve models,'' Stat. Med. 35, 2031--2057 (2016).
//
inline double AUC_binormal_LR(
                              const double d_a,
                              const double c
                              )
{
    // NOTE: The BVN part of the calculation SEEMS to be only a (small) correction to the CDF of the normal distribution.

    const double inf = 10.; // large number, taken to be infinity for the BVN integrals
    const int    n   = 64;  // number of Gauss--Legendre points

    //=========================================================

    double AUC;

    // CALCULATE CDF OF STANDARDIZED BIVARIATE NORMAL DISTRIBUTION
    //
    // NOTE: There is not a closed form expression for this.
    //
    // NOTE: ... See:
    //
    //     https://stats.stackexchange.com/questions/48227/correlated-bivariate-normal-distribution-finding-percentage-of-of-data-which-is
    //
    // NOTE: This is therefore calculated by numerical integration.

    double rho = -(1. - c*c)/(1. + c*c);

    double b1 = -d_a/std::sqrt(2.);
    double b2 = 0.;

    std::function<double(const double)> f_outer = std::bind(
                                                            CDF_standardized_bivariate_normal_distribution_outer,
                                                            // -----
                                                            std::placeholders::_1,
                                                            // -----
                                                            std::cref(rho),
                                                            // -----
                                                            std::cref(inf),
                                                            std::cref(b2),
                                                            // -----
                                                            std::cref(n)
                                                            );

    double F_BVN = functionsxx::integration::quadrature_GaussLegendre(
                                                                      f_outer,
                                                                      // -----
                                                                      (-inf),
                                                                      b1,
                                                                      // -----
                                                                      n
                                                                      );

    F_BVN *= ( 1./(2.*PI*std::sqrt(1. - rho*rho)) );

    // CALCULATE AUC

    // NOTE: See: Metz and Pan, or Hillis Eq. (27).
    AUC = boost::math::cdf( boost::math::normal(), (d_a/std::sqrt(2.)) ) + 2.*F_BVN;

    // RETURN
    return AUC;
}


//
// DESC: Outer integral for calculating CDF of standardized bivariate normal distribution.
//
static double CDF_standardized_bivariate_normal_distribution_outer(
                                                                   const double x,
                                                                   // -----
                                                                   const double rho,
                                                                   // -----
                                                                   const double inf,
                                                                   const double b2,
                                                                   // -----
                                                                   const int    n
                                                                   )
{
    // DO INNER INTEGRAL (FOR FIXED x)
    std::function<double(const double)> f_inner = std::bind(
                                                            CDF_standardized_bivariate_normal_distribution,
                                                            // -----
                                                            std::cref(x),
                                                            std::placeholders::_1,
                                                            // -----
                                                            std::cref(rho)
                                                            );

    double F = functionsxx::integration::quadrature_GaussLegendre(
                                                                  f_inner,
                                                                  // -----
                                                                  (-inf),
                                                                  b2,
                                                                  // -----
                                                                  n
                                                                  );

    // RETURN
    return F;
}


//
// DESC: (Integrand) of CDF of standardized bivariate normal distribution.
//
static double CDF_standardized_bivariate_normal_distribution(
                                                             const double x,
                                                             const double y,
                                                             // -----
                                                             const double rho
                                                             )
{
    return ( std::exp(-(x*x - 2.*rho*x*y + y*y)/(2.*(1. - rho*rho))) );
}
